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ABSTRACT 

A control volume method is proposed for planar div-curl systems. The method is inde- 
pendent of potential and least squares formulations, and works directly with the div-curl 
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and two control volumes rather than the other way round. A discrete vector field theory 
comes quite naturally from this idea and is developed in the paper. Error estimates are 
proved for the method, and other ramifications investigated. 
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1 Introduction. Although div-curl systems occur in fluid dynamics, in elec- 
trodynamics and several other applications, relatively few discretizations are 
available. A possible reason is that equivalent potential formulations often ex- 
ist. Another reason is that the simple minded Galerkin finite element approach 
is not convergent in general. A least squares formulation is the usual way to 
get a convergent finite element scheme for planar problems. The situation in 
three dimensions is worse. For example, the vector potential approach can have 
spurious mode problems [10], Underlying the difficulties is the latent overdeter- 
mination of the div-curl system (4 equations, 3 unknown functions). 

This article contains an alternative approach which works directly with the div- 
curl equations and does not involve potentials or least squares. The new recipe 
is finite volume based, but differs in a key ingredient. The change permits 
the development of a discrete vector analysis. In turn, this provides tools for 
analysis of the discretization and for other purposes. 

Central to the approach is the use of dual pairs of meshes made up of “comple- 
mentary volumes”. In three dimensions, they have the property that the edges 
of each mesh are perpendicular to the faces of the other. There are many possi- 
bilities for such mesh pairs. In two dimensions, the simplest example consists of 
the staggered Cartesian meshes (MAC meshes) well known in fluid mechanics. 
The complementary volumes are the basic mesh squares and the shifted mesh 
squares centered on the nodes of the basic mesh. For triangular and tetrahedral 
meshes, an example is given by Voronoi- Delaunay mesh pairs. Many other pos- 
sibilities exist, including prismatic meshes in three dimensions .and combinations 
of these meshes. 

The basic idea of the discretization is to define field components along the edges 
of one of the meshes, and, therefore, normal to the faces of the_othe_r. In two 
dimensions, this single component is enough to permit the definition of two 
field operators corresponding to div and curl. With boundary conditions, these 
are sufficient to define the discrete field. Associated with the nodes of the two 
meshes are the two discrete potentials which generate the null spaces of the 
discrete div and curl. The usual relations are valid between these operators and 
potentials. 

We will address only the two dimensional problem in this paper. The main ideas 
go over naturally to three dimensions but are sufficiently different to warrant a 
separate treatment. It will be given in a forthcoming report. The goals here are 
to provide a framework for error estimation of complementary volume schemes, 
and to develop the main tools of the discrete vector field theory. 

In sections 2 and 3 the class of meshes of interest is specified, and a typical 
discretization is derived. Sections 4 and 5 contain the formulation and prop- 
erties of the discrete vector field theory and some applications. Sections 6-8 
are concerned with the main results of the error analysis. Section 7 in partic- 
ular establishes a link with finite elements and applies the previous results to 
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error estimate a potential theory problem. Section 9 discusses special topics 
concerning rectangular meshes. Section 10 extends the results to more general 
boundary conditions. 

The discretization technique reported here has significant extensions to higher 
order systems of partial differential equations, particularly to viscous fluid flow 
problems. Algorithms for the Navier-Stokes equations are provided in [4] and 

[7]- 


2 Locally Equiangular Meshes. We begin by defining a class of meshes of 
interest. Let Cl denote a bounded polygonal region of R 2 . Cl may be multiply 
connected. Let Tq denote the outer boundary, let i = 1,2,. . . , r denote the 
inner (polygonal) boundaries if they exist and let T := Uj_ 0 IV T* i — 1, 2, , r 
are considered to bound “holes” Q* t = 1, 2, . . . , r . Cl will be triangulated and 
a dual tesselation will also be used. Let r denote a triangulation of Cl with T 
triangles denoted by T{ i = 1,2, .,T with N nodes X{ i = f , 2, . . . , N\ € Cl 
and xj j = Ni + 1 , N\ + 2, . . . , N E T as well as E edges cq i = 1 , 2, . . . , E\ 
E Q ( interior edges) and <jj j = E\ + 1, E\ -{- 2, . . . , E E T ( boundary edges). 
Two triangles are called adjacent if they share a common side. Then we can 
construct a dual tesselation by joining the circumcenters of adjacent triangles. 
Many duals of a given triangulation can be constructed. For example, another 
one could be made by joining the centroids of adjacent triangles. The dual which 
is based on circumcenters is rather special and will be called the normal dual 
since if two triangles are adjacent, the line joining their circumcenters is normal 
to (and bisects) their common side. The dual figures are polygons and in general 
they can have self intersecting boundaries. They will be called covolumes . The 
covolume associated with an interior node is the polygonal figure obtained by 
joining (in order) the circumcenters of the adja-cent triangles w_hich share it. 

We will associate a (boundary) covolume with each boundary node as well. The 
procedure is illustrated in Figure 1, where the covolume for the boundary node 
A is the interior of the polygon PATSRQ. 

In Figure 1 (3,R, and S are the circumcenters oF their triangles, and P and T 
are the midpoints (and circumcenters) of their edges which are on F . 

The normal dual, consisting of T nodes, E edges and N covolumes is denoted 
by r / . Sometimes, we will use the “co” prefix to denote various elements of the 
normal dual, for example in referring to coedges, comesh and so on. 

There is additional complexity associated with self intersecting covolumes which 
we wish to avoid. To do so, we require that r is “locally equiangular” [8]: 

Definition. A triangulation is locally equiangular iff for every pair of adjacent 
triangles which form a convex quadrilateral, the sum of the angles opposite the 
common side is at most 180 deg. 


2 



Elementary geometry shows that if r is locally equiangular then (1) each interior 
covolume is convex and (2) distinct covolumes have empty intersection and each 
point of Q is either in a covolume or on the common boundary of two covolumes. 
Then r' has inner and outer coboundaries V consisting of those coedges which 
intersect edges of r with just one node on T. These will be denoted by r' 0 , and 
r' i = 1, 2, . . . , r when the latter exist. The section QRS in Figure 1 is part of 
T*. Note that any part of T' might penetrate the corresponding part of r . This 
will happen if there are obtuse angles opposite an edge of T , even if r is locally 
equiangular. 

We will obtain the results assuming that r > 0. The modifications for r = 0 are 
mostly self evident. 

The two tesselations r and r' are close to being a Delaunay- Voronoi pair. How- 
ever, although a Delaunay triangulation is always locally equiangular, the con- 
verse is false, at least when the classical definition of the Voronoi diagram is 
used. The problem occurs at boundaries. Recall that a standard Delaunay tri- 
angulation is defined by joining adjacent vertices if their Voronoi figures share 
a common edge. Then it follows that the Delaunay triangulation is a trian- 
gulation of the convex hull of the vertices. Since Q is not convex in general, 
constructing the Delaunay triangulation of the vertices of r does not necessarily 
give a triangulation of ft . For the purposes of this article these distinctions are 
unimportant. Local equiangularity is the only property we need since all the 
properties required can be derived from it. 

Some of the results below must be modified for certain kinds of trivial triangula- 
tions, in particular those with no interior triangles. We will always assume that 
the triangulations are sufficiently fine for the purposes at hand. No significant 
loss of generality is incurred by this assumption. 

We will make frequent use of the following special case of Euler’s formula: 
Lemma 2.1. Let r denote a plane triangulation with N vertices, T triangles, 
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E edges, and r holes. Then 


N+T— E+l — r. 

Proof. If r = 0 the lemma is easily proved by deleting triangles from r while 
maintaining a count of N, T and E. If r > 0 we can imagine the holes trian- 
gulated consistently with r . This gives a triangulation say t without holes. 
Combining the already established result for f and for the triangulations of the 
holes, the result follows by subtraction. □ 


3 Div-Curl Systems. One div-curl system of interest is 


div u = p 

in Q 

a) 

curl u = w 

in 

(2) 

u n = / 

on r 

(3) 

= 7i i = 


(4) 


/ U '' 

J i\ 

where u .= (u, u), curl u := v x — u v , t denotes the positively oriented unit 
tangent, n the outward unit normal, and it is assumed that 


J pdxdy = J fds. 


(5) 


Also assumed is that p € L 2 (Q), u e L 2 (fi) and / 6 ff»/»(r) and that the 
system has a unique solution u 6 See [6], [9] for information on this 

point. We will explain the basic ideas of the discretization in terms of (l)- (4). 
Section 10 extends the results to other boundary conditions. 

Referring to Figure 2, we approximate (1) by 


Uihi + «2^2 + U3/13 = / pdxdy 
in 


( 6 ) 


Here and below, the u, denote approximations tou nj, n ; denote unit normals 
directed outwards and hj > 0 denote the ordinary side lengths of r. There will 
be a similar equation for each one of the T triangles in r . In matrix form, with 
u denoting the vector of components u , , these flux equations can be written as 

Fu = p (7) 

where u E R E and p £ R T . Analogous to (6) discrete fluxes from the holes are 
also defined and we denote them by 


n«;r>) 


i = 12, ..., r. 
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For any u E R E , it is also convenient to introduce the r x E matrix T and 
denote the hole fluxes by Tu . 

To approximate (2) we integrate it over an interior covolume r[ E t* as illus- 
trated in Figure 3. The arrows on the covolume edges denote the directions of 



the unit normals to the associated triangle edges. Approximation of the integral 


/ u • t ds (= curl u dx dy) 

Jdr[ Jri 

where denotes the positively oriented boundary and t the unit tangent gives 

Y^Ukh'u = / w dx dy (8) 

l vrj 


where the sum is over the coedges of and h'^ > 0 denotes the length of a 
coedge. Assembling these circulation equations gives the matrix equation 

C 0 u = w ( 9 ) 
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where u £ R and u £ R N ‘ . The zero subscript is intended to suggest that 
circulations in (9) are computed for interior nodes only. Discrete circulations 
around the hole coboundaries are defined similarly to (8) and denoted by 

j = 1,2, 

These circulations will be represented as Cu where u £ R E and C is r x E. The 
boundary condition (3) is discretized by defining boundary edge values for u by 

Uk = r k l fds k = Ei + (io) 

Note that there are N — N x of these equations. 

The prescribed circulations (4) are approximated as 

C( u ;r>) = Tj + I u> dx dy j — 1,2, ...,r (11) 

Jsj 

9j j = 1, 2, r 

where Sj denotes the strip lying between Tj and T'- . We will assume that T does 
not penetrate to avoid extending u outside fi . This will hold iff opposite 
every triangle edge on Tj j — 1,2, ...,r is an acute angle. (No restriction is 
necessary for a simply connected domain). There is, of course, no assurance that 
a locally equiangular or even a Delaunay triangulation has no obtuse angles. 

Equations (7), (9), (10) and (11) are a linear system of T+ Ni + (N - N^ + r = 
T + N + r equations in E unknowns. By Euler’s formula there is one more 
equation than unknowns, so we expect a single constraint on the data This will 
turn out to be 

^2 * = J2 A ( 12 ) 

<7 je r 

which holds by (5). We will prove in section 5 that these discrete equations 
indeed have a unique solution. 


4 Mesh Matrices. In this section some basic properties of dual mesh systems 
are derived. These will lead in the next section to a more detailed formulation 
of the equations m section 3. Similar results are valid for more general mesh 
systems than plane triangulations but will not be needed below. 

Let r denote an arbitrary triangulation of with T triangles, E edges including 
E x interior edges and N nodes of which N x are interior nodes. Label the interior 
nodes 1,2,.. N u the interior edges 1 ,2,...,£’ 1 and assign the positive direc- 
lon along each edge <r k k = 1 2, . . . , E to be from lower to higher node number. 
Let t denote a dual mesh, with T nodes E edges and N covolumes. The dual 
can be quite general, subject to having exactly one node in each triangle. The 
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dual edges obtained by joining adjacent dual nodes are in a biunique correspo- 
nence with the edges of r. For each boundary triangle, we join the dual node 
to a point on the triangle’s boundary. The dual edge corresponding to cr* is 
denoted by cr' k . For orienting the edges of r’ we use the convention that (cr-^cr*) 
are oriented like a Cartesian coordinate system in the plane. This convention 
applies to both interior and boundary edges. 

Denote by D the E x N edge-node incidence matrix of r, where 

{ +1 if <T{ is directed into node j 
— 1 if cr j is directed out of node j 
0 does not meet node j . 

D is easily seen to have rank N — 1. Define Do to be the E\ x N\ matrix 
obtained by deleting rows and columns of D corresponding to boundary edges 
and boundary nodes respectively. Do has rank N\. If r > 0 we also define an 


Ei x r matrix V as follows. 


r +i 

if < 7 ; is directed into T s 

Vi, := { -1 

if a i is directed out of T s 

1 o 

a i does not meet T a . 


In this, cr, denotes interior edges, and s = 1,2, . . . ,r. V has rank r. 
For the dual r f we define the E x T incidence matrix B as follows: 


Bij 


-hi if o\ is directed out of Tj 
— 1 if is directed into Tj 
0 (J x does not meet Tj. 


B has rank T. B 0 denotes B with rows corresponding to boundary edges deleted. 
Bq is of order E\ x T. If r > 0 we define an E x r matrix B by 


Bu : = 


-fl if a[ is directed out of 
-1 if (7- is directed into 
0 <j,- does not meet T a . 


B has rank r. 

An important result is the following: 

Theorem 4.1. Let v E R El . Then 3<f> E R T such that v = B Q <f) iff D f 0 v = 0 
and V*v — 0. 

Proof. Since rank (Do) = T — 1, we have 

dimA r (Dg) — 
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Ei-T+1 

N x +r 



using lemma 2.1. Here and below, N(-) and R( ) denote the null space and range 
of their arguments. Next, it can be verified directly that BqDqiP = 0V^> G R Nl , 
and BqVZ = 0V£ € so that i?(D 0 ) C jV(Bj) and i?(P) C N(B J). Direct 
verification also shows that i?(Do) PI i?(D) = 0. Since dim R(Dq) -f dim i?(2>) = 
Ah + r it follows that JV^f^) = i?(D 0 ) U #(£>). Solvability of the equation 

Bq<^> - v 

holds iff 

(v,z) = 0 Vza®, 

( , ) denoting the standard Euclidean inner product. By the above, this is 
equivalent to 

(■ v,D 0 ip ) = 0 vV'eft" 1 

Mo = o VUR T , 

and these in turn are equivalent to the theorem. □ 

Corresponding to theorem 4.1 is: 

Theorem 4.2. Let v £ R E . Then 3^6 R n such that v = Drp iff B*v = 0 and 
B l v = 0. 

Proof. Similar to theorem 4.1. □ 

5 Discrete vector fields. In this section we will express the flux and circu- 
lation operators in terms of the mesh properties of section 4 and develop some 
analogs of vector field theorems. For the remainder of the paper we will be 
using the normal dual exclusively. In any normal dual, the distance between 
two circumcenters will be zero if they coincide. Although this situation is “non- 
generic” , it does occur in some situations (section 9). Other than this, it is 
assumed throughout that the circumcenters are all distinct. If, in fact, there 
are coincident circumcenters, the results obtained remain true, but minor varia- 
tions in some proofs may be required. A second assumption, made throughout, 
is that the coboundary F' does not penetrate T. This was already mentioned in 
section 3 following (11). 

In R e introduce the inner product [•, ■] defined by 

[u, v] := (u, HH'v) = (u, H'Hv) (13) 

and denote the associated norm by 

II « lk= [«. «] 1/2 - 

In (13) H := diag(/ifc) and H ' := diag(/i] b ). Both H and H f are invertible and 
we let W H H ( . In Figure 4, ABC and ACD are adjacent triangles from r, 


8 


in niMiaiiilfillKlIii iiii i Imhiiii in| I Hill 



A 



Figure 4: 


and P and Q are their circumcenters. Let | AC |:= /i* and | PQ |:= h k . The 
area of the kite shaped figure APCQ is hkh! k /2. The corresponding areas at 
boundaries are /ij/iJ/4. It follows that || • \\w is twice a discrete L 2 (Q) norm. 
This interpretation holds for any locally equiangular triangulation. 

Denote R E equipped with [•, •] by U = U(Cl). Then we can refer to “grid 
functions” u E U(Q) and regard them as having boundary values u[r and 
interior values u\n. Define 

Uo := {ue U]u\ r = 0}. 

The flux and circulation operators of section 3 can be expressed as 

F = B'H 
Co = DlH'o 
C = V l H' 0 
T = B t H 

where H$ denotes the restriction of H f to interior edges.. Verification of these 
is by direct calculation. Note that orientations are automatically taken into 
account in this formulation. Difference operators are defined as follows: 


s<t> 

CQ 

1 

11 

v<£e r t 

So<f> 

:= (H' 0 )~ 1 Bo<f> 

V<£G R t 

Txl> 

:= H~ 1 Dip 

W< G R". 


We will the notation Ho to denote the restriction of H to interior edges, and 
Wo := H 0 H' 0 ~ #o#o- The theorems of section 4 translate into theorems about 
the existence of potentials for mesh functions u with Cqu = 0 and Cu = 0 
(velocity /scalar potential) and Fu = 0 and Tu = 0 (stream function/vector 
potential). 

Theorem 5.1. If Cqu = 0 and Cu = 0 then there exists <j> E R T such that 
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u = So<j). Conversely , if u So<p then Cqu = 0 and Cu = 0. 

Proof. Use theorem 4.1 to prove the first part and direct substitution for the 
converse. □ 

For the stream function, the analogous result is: 

Theorem 5.2. If Fu = 0 and Fu = 0 then there exists rf E R N such that 
u = T / tp. Conversely, if u := Tip then Fu = 0 and Fu = 0. 

Proof, For the first part use theorem 4.2. The converse is easily proved by 
direct substitution. □ 

Note that the converses in these theorems furnish analogs to the vector identities 
“curl grad <p = 0” and “div curl u = 0” . 

Another useful result is: 

Lemma 5.1. For all u E U , <j> E and ip E R N we have 

(1) [u, S<f>] = (Fu, <P) 

( 2 ) [u,m = {Cu,i>). 


Proof. (1). By definition of S 

[ti, S<P ] = (u, HH\H ( )- l B<P) 

= (B t Hu,<P) = (Fu ) <P). 


The proof of (2) is similar. □ 

These are analogs of integration formulas. For example, the first is analogous 
to 

/ u. V <t>dx dy = — 

Ja 

Another useful result following from above is an analog of the Helmholtz de- 
composition of vector fields. In general, this decomposition is for the subspace 
of Uq C Uq whose definition is: 

U 0 := {wE Uo;Cu = 0 and Fu — 0} 

Actually, since u E Uq ^ Fu = 0, this requirement in the definition is redun- 
dant. It is included to make explicit a symmetry in the hypotheses. Also, before 
computing Cu we must restrict u to interior edges. This point will arise several 
times below and also in connection with Cq. We will still denote the restriction 
by u , and without explicit mention each time. If r = 0, then Uq = Uq and the 


/ <pdivudxdy (<^> |r= 0). 

Jrt 
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decomposition is for U 0 itself. We also define the following subspaces of U 0 - 

Z 0 = {« € Uq\ Fu -- 0} 

Wo = {«€ U Q \ C 0 u = 0}. 


Now we have: 

Theorem 5.3. Uo has the decomposition 

Uq = Zo © Wo , 


which is orthogonal relative to [*, •]. 

Proof. First, note that if w € Wo, then by theorem 5.1 and part (1) of lemma 
5.1 Vz € Uo, 

[z, w) = [z, So<i>] = [z, S<j>] = (. Fz , 4) V<t> e F T - 

Hence, if z is orthogonal to Wo, then Fz — 0 which implies z € Zo- O n the 
other hand, if some z € Uo satisfies Fz = 0 and also C 0 z = 0 then we have for 
some 4>' G R T 


[z, z] = [z, S<t>'] 

= (Fz,4') 

= 0 


so that z = 0. This proves the result. □ 

As an application of theorem 5.4 we will prove that the discrete system of section 
3 has a unique solution. These equations are 


Fu = B*Hu = p 

(14) 

C 0 u = DlH'oU = w 

(15) 

Cu = VH' 0 u = g 

(16) 

Xu = /. 

(17) 


Here, u G R E , p, u, f and g are as before and the E x (N - Ni) matrix I 
denotes the identity restricted to the boundary edge values. 

Theorem 5.4. Equations (U)-(l 7) have a unique solution u G U. 

Proof. Consider the homogeneous problem Fu = 0 Cqu = 0 Cu = 0 lu = 0. 
In this case, u G W 0 f\ Z 0 and by theorem 5.3 it follows that u = 0, giving 
uniqueness. Consequently, the matrix [F 1 Cq X 1 C‘] which is of order (E+ 1) x F 
has rank E. The augmented matrix 

' f* c* 0 r C % 1‘ 

. ? & P ? \ 
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of order (E + 1) x (E + 1) has rank at most E. This follows by subtraction of 
the third block row from the sum of the first block of rows followed by use of 
the compatibility condition (3.12). Hence, existence follows. □. 

A method for solving (14)-(17) by reducing them to potential equations can 
also be found using the results above including theorem 5.3. We can suppose 
without loss of generality that / = 0, simply by substituting the known values 
of u into (14) and (15). These values do not appear in (16). We can also assume 
g — 0 since it is trivial to find a solution of (16) and then modify u accordingly. 
Hence, it is only necessary to consider (14)-(17) with / = 0 and g = 0. That 

means u £ U 0 . Now seek u as a sum «('> + t/*), where both of these are in U 0 
and 


and 


Ftfl> = p 
Con® = 0 


= 0 

C 0 u W = Q. 


(18) 


Consider the first set of equations. By theorem 5.1 we can write u^) = an( j 
(18) then becomes f 

FS 0 <j> = p. 

Taking account of the boundary conditions, this can be factored into 

[5‘ W o 1/2 ][W o 1/2 5 o ]0 = P- 

The coefficient matrix is positive semidefinite since B 0 has rank T- 1. Clearly, 
this reduction parallels a procedure for the continuous problem. The second 
set of equations can be solved by a similar approach. The details are given 
(for a different context) in section 7. These procedures reduce the equations 
(14)-(17) to the solution of positive semidefinite equations. Various iterative 
procedures for the solution of (14)-(17) turn out to be disguised iterations for 
these equations. 


6 Error analysis. In this section we will estimate the error in approximating 
the solution u of the div-curl system (l)-(4) by the solution u of the covolume 
approximation (14)-(17). 

To begin, we need the following result: 

Lemma 6.1. Assume that v £ Wo and that u £ Uq and Fu = 0 Then 

[«, v] = 0. ' ' ' ‘ 

Proof. In this lemma we do not require that u £ U 0 . By theorem 4.1, and using 
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u| r = 0 we have [«, v] = [u, So*] = [«, S*] = (Fu, *) = 0 using lemma 5.1(1). 

□ 

Next, we introduce the following “mesh functions” u (,) G U i = 1,2 computed 
from the exact solution u of the div-curl system as follows: 

:= — f u n ds k = l,2,...,E (19) 

hk Ja k 

where cr k is traversed positively, and n points to the right of cr*, and 
u™ := tt f n-nds fc = 1,2, . . . ,Fi 

h k M 

u™ := k = E^ + l,...,E 

where cr' k is traversed positively, and n points along a' k . 

Now denoting £^ := u — v£'\ i — 1,2 we have 

F(^ = 0 eW e U 0 

<V 2) = o £< 2 > G Uo 

Then by lemma 6.1 

[ € (D )f ( 2 )] = 0 

and consequently, defining u := (u^ + t/ 2 ^)/ 2, it follows that 

[„-a, u -5] = [i(.< 1 > + « <,) ). !(»<'> + < ra )l 

= I[ C (D _ e ( 2 ), f (D _ t (2>] 

= IfuC 1 ) - u (2 >, u (1) - U (2) ]- 
4 l 

The right side is independent of u, while depending only on u so this gives a 
preliminary evaluation of the error. We can estimate the left side in terms of 
for example, by means of 

|j u — u ||w = ||„-«W-( u ( a )-«W)/2|k 

> || u-u (l) llw II u (2) -u (1) \\w 

so that the estimate becomes 

II u - u (1) ||w < II « (2) - « (1) Ik • ( 2 °) 

Evaluation of the right side of this proceeds as follows: referring to Figure 5 let 
I< denote a kite shaped domain with perpendicular diagonals s and s. 
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Figure 5: 

We shall assume the diagonals are aligned with a Cartesian coordinate system 
intersecting at the origin 0. Let the diagonal s on the i-axis extend from x = —L 
to x — L and let the diagonal s' on the y-axis extend from y = — Mj to y = M 2 . 
Define linear functionals //’) i = 1,2 on by 

1 

^ (1) ( v ) := 2L J L V2 ( X ’ dx 

1 

" (v) := 

where v := (v lt v 2 ) 6 H By the trace theorem, //(*)(.) j = i 2 are 
bounded onH 1 (K) and so is - // 2 >Q. Also, if v is constant in K, then 

i/i /i )(v) _ o. It follows that there exists a constant C(K) such that 

| (M (1) - ^ (2) )(v) |< C(K) I V \ 1>K 
where [ • ^ denotes the seminorm. 

A standard scale change argument shows that C(I<) depends on the ratio 


max 


(- -) 

\ L 1 M ) 


where M | M x + M 2 |. In terms of mesh geometry, for the kite associated 
with crjc this becomes 

fh k h'\ 

Substitution into the right side of the error estimate now gives 

,,m - Ilk = £ I I // 11 -,, n >)(v) | ! | h,h' t I 


< Ecw i,A-*max(*t, (^i) 2 ) 

k 

< C ma x(h 2 , h' 2 ) | u \\ n 
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where h := max^ hk and h f := max*; h f k . This proves: 

Theorem 6.1. Assume that u E H 1 (Q). Then with u denoting the approximate 
solution and u ^ computed from (19), we have the estimate 

|| u — || w < Cmax(/i, h') | u |i f n ♦ 

Note that there is implicit dependence on the angles of the triangulation through 
the appearance of the dual mesh parameter h* . 

This estimate can be improved for certain regular meshes. The key observa- 
tion is that the functional ~ vanishes on linear functions in K 

if the vertical diagonal in Figure 5 is also bisected. This will occur in a gen- 
eral triangulation iff all the triangles have the same circumradius. A complete 
characterization of such triangulations is not known. It certainly includes the 
standard uniform triangulation of the unit square, and the more symmetric 
triangulation in which each mesh square of the uniform rectangular mesh is 
subdivided by its diagonals. It also includes triangulations made of equilateral 
triangles. (The first two of these examples are brought within the covolume 
framework in section 10.) To exploit the mesh regularity, we assume now that 
v E H 2 {K). Then it follows that 

K/iW-^Kv) |< C(I<) I V \ 2 ,K 

and the constant depends on 


< L3 

max( — , 
M 


M 3 
L ’ 


ML). 


Proceeding as above, we obtain for the regular meshes the estimate 
|| u - u (1) ||tv < Cmax(h 2 , (ti) 2 ) \ u | 2 ,n • 


From the finite difference viewpoint, this means that the covolume scheme is 
second order accurate. The analysis strongly suggests that rapid changes in the 
circumradii of adjacent triangles should be avoided in practice. Similarly, we 
expect that meshes which vary smoothly in this sense will yield better accuracy. 
It would be possible to estimate the deviation from second order accuracy in 
terms of the change in adjacent circumradii, but we omit this. 


7 Covolume solution of Poisson’s equation. An interesting application 
of the stream function can be given. Consider the equation 


-A* 

* |r 


e H 2 (n)nH^(Q) 


( 21 ) 

( 22 ) 


0 . 
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Integration over a covolume rj gives 




and using the approximation Tip for we get after approximation of the line 
integral 

—CqT ip = ah 

In more familiar notation, referring to Figure 6, a typical component of this 
equation is 

(* 3 > 

Boundary conditions are used in the obvious way. We assume a locally equian- 



h'i 


Figure 6: 

gular triangulation. On rectangular meshes, and to a lesser extent on triangular 
ones, this approximation is well known. 

Assembling the equations similar to (23), we obtain a linear system 

Kxp = Q 

where K is of order N\ x N\. K has the unexpected property of being the 
same as the piecewise linear finite element matrix for (21)-(22) on r. That is, 
if Xi i — 1, 2, . , , , N\ denote the standard piecewise linear basis functions for r, 
then 

K ij- (\7\i)(v\j)dxdy i, j = 1,2, . . . , Ni. 

Jn 

A proof may be found in [1]. Another proof could be based on the potentially 
useful factoring of K which follows by letting G := Hq 1 Do. In terms of this we 
have 

K = GHVqG. 
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This shows immediately, for example, that K is positive definite since Dq has 
rank N\. 

A difference from the finite element formulation is apparent in the formation of 
the data vector Q. The components of this vector are formed by integration of u > 
against the characteristic functions of the covolumes. This is in contrast to the 
integration against the A , called for in exact finite element theory. It is possible 
that an error estimate for the covolume scheme could be derived from the finite 
element estimate by the usual techniques for handling quadrature formulas. On 
the other hand, it seems quite natural to obtain an estimate within the covolume 
framework. This can be done as follows. Define it := Tty . Then we have 


Fu 

= 0 


CqU 

= w 


u | r 

= 0 


Cu 

= 7 i 

i - 1, 2, r 


where 7,- is computed as indicated. The exact solution of the analogous div-curl 
system is u = (t/i , 1^2) where U\ = d y ty and it2 — — d x ty. Based on u define u ^ 
by (19). Since u | r = 0 and Fu W = 0 by theorem 4.2 there exists E R N 
such that = Tip^\ Next, we have 

Lemma 7.1. 

I <=*(*<) i= 1,2,... , AT. 

Proof. By definition, integrating positively along <7*, we have for k = 1, 2 , . . . , E 
iti 1 ^ := 7- / u nrfs 

hk Jau 

= Y j - d *V) nds 

- Us* 

= (T*)| fc 

so that |jb= (T'tp^) [fc. Recalling that rank(T)=rank(D)=A r — 1, it follows 

that ty(xk) and ip(x *) differ by a constant which can take to be zero. □ 

Define 

a($, tf) := [ (v^)(V^) dx dy 6 

The main result is: 
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Theorem 7.1. Let xp denote the piecewise linear interpolant of xp the covolume 
approximation to 'F. Then 

|| - V ||i,n < C max(/i, ti) || V || 2 ,n • 


Proof. By theorem 6.1, we can write 

||u-u (1) |k = || Tip - Tip (1) ||jy 

< Cmax(ft, h') || ^ || 2 ,n 

and so 

(C 0 T(ip - ip w ), ip - ^ (1) ) < C[max(/i, h ') || $ || 2> n] 2 - 
By the lemma, denoting by ^ the piecewise linear interpolant of on r, 

a(xp - xp — #) < C[max(ft, h f ) || || 2 ,n] 2 

and 

a(xp - , $ - tf) < C[max(ft, ft') || || 2 ,n] 2 + a(# - - tf). 

Using approximation theory it follows that 

a(i/> - $ y xp — $) < C[max(ft, ft') || || 2 ,n] 2 

and the result follows from Poincare-Friedrichs’ inequality. □ 

Thus, the covolume scheme and the finite element scheme are of the same order 
in the H l (Q) norm. 


8 Tangential Components. The discretizations of section 3 produce approx- 
imations to the velocity components normal to the triangle edges. Tangential 
field components are known in the directions of the comesh edges. This is a 
basic characteristic of the covolume method. But in many applications it is a 
vector field that is required and not merely sets of components. An example 
occurs in connection with approximating convection terms in viscous flow prob- 
lems [7]. In this section we give an algorithm for computing a set of tangential 
components from a given set of normal components. In this way, a (discrete) 
vector field is obtained. The covolume scheme itself is not changed. Instead, the 
tangential components are found from the covolume solution in a “postprocess- 
ing” operation. Analysis of the error of the tangential approximations is also 
presented in this section. 

We begin by noting that the three normal field components for a given triangle 
are too many to determine a constant field in the triangle. The following result 
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gives the compatibility condition for the three components to uniquely determine 
a constant field. 

Lemma 8.1. Given three normal components i = 1,2,3 on triangle sides 
with lengths hi and corresponding unit outward normals n» where i = 1,2,3 
there exists a constant vector u such that 

u * n; = Vi i = 1,2,3 (24) 


iff the flux 


V\h\ + t>2 ^2 d~ — 0’ 


u is unique . 

Proof. Equation (24) is a system of three equations for two components of u. 
For nondegenerate triangles r the coefficient matrix A (ni)j clearly has rank 
2, so it follows that dim# (A*) = 1. Now the equation 

Aw = 0 


certainly has the solution w = (hi, hi, h 3 ), since 

0 = jv(l)dzdy = L n ds = nj/t, + n 2 h 2 + 113 / 13 . 

Hence, (24) will be solvable uniquely iff 

V\ hi + v 2 h 2 + V3/13 = 0 
which is what we wanted to show. □ 

From this lemma it follows that if Fu = 0, then we can construct a piecewise 
constant vector field u on r whose (continuous) normal components are given 
by u. An easy way to compute a tangential component along an edge of the 
triangulation is to average the tangential components of the constant fields in 
the triangles sharing the edge. The remaining problem is then to construct 
the piecewise constant field when the flux is nonzero and lemma 8.1 does not 
apply. There is no unique way to do this. For example, we can take any pair 
of normal components and designate them as the respective components of the 
sought vector in a triangle. This procedure will be used. Denoting the kite areas 
associated with the triangle by A'i > I< 2 > A 3 , we choose the corresponding 
components u, and u 2 associated with the largest and second largest kite areas 

to define the vector. 

To compute a constant interpolating field u/, let u := («i, «p) and define 
cos# := ni • n 2 - Representing u / := ouni +or 2 i 1 2 gives the equations 

( 1 cos 9 \ f 0-1 \ _ f «1 \ (25) 

^ COS 9 1 ) \ <*2 J \ u 2 J 
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which uniquely determine oq, (* 2 ■ It can be checked from this representation 
that 

^ 2 „ _ 

u i u i < — 27 « * ( 26 ) 

sin 9 v 1 

a result we will use below. 

Interpolation of the u/’s across the common edge of their triangles is also subject 
to a certain amount of arbitrariness. One possibility is to use the vector from 
one of the triangles and ignore the other. Another is to weight the vectors 
equally, while a third is to use full linear interpolation. In this last case, it is 
necessary to assign a localization of the vectors. The most natural choice is 
to attach each u / to the circumcenter of its triangle. To justify this, we can 
observe that the two directions from which u/ is made indeed pass through the 
circumcenter. The line joining the circumcenters is divided in a certain ratio by 
the common edge of the triangles which determines the interpolation weights in 
the usual way. Note that if one of the triangles is obtuse, this becomes a linear 
extrapolation^ Each of these methods corresponds to a weighting of the form 
Aiu 7 -f A 2 U 7 in which A x + A 2 = 1. Some situations suggest other weightings 
to use. At boundaries in particular, a (1, 0) weighting is required. In other cases 
the weights can vary from edge to edge. 

Next, we will estimate the error of this tangential approximation scheme. The 
interpolated vector for the k th edge has the form 

wW:=Af , ^ l) + ^ ) u f 3) 

where the superscripts (&i), (£ 2 ) denote the two triangles on each side of the 
k f edge. The tangential component of this is just For u E U, let 

Tu E U denote the mapping to the tangential components generated in this 
way. 

A parameter which will appear in the estimate is the ratio \J I<[^ f =: ^ k \ 

where are the kite areas for tj. We will define S := ma This 
quantity will be estimated later. Also we will define 0 1 = min m (|sin $( m )|) 
where 6 ( ^ corresponds to the angle in (25) above and the minimum is over the 
triangles. In addition to these we will denote by A the maximum absolute inter- 
polation coefficient: A := max^A^], |A^|). The next result will be required 
below: 

Lemma 8.1. The following bound on T holds: 

\\Tu\\ w < Vu€C. 

Proof. Let denote the kite area on the edge c ^ . Then since 

| t W . w (*)|2 < 2 A 2 (|| 4 *' ) ||?, + || uJ *» ) || i 3 ,) 
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we have, using (26) 

k k 

k 

Also 

|«(*‘)| 2 A<*> = [(4 10 ) 2 + (u ( 2 fcl) ) 2 ]A’ (i) 

< 6 2 [( U ( 1 l ‘ ) ) 2 /^ fc,) + (4 i ‘ ) ) 2 4 t,) ] 

with a similar bound for |«(^)| 2 A'( |s) . Noting that each product appears at most 
four times and substituting these into the previous inequality gives the result. 
□ 

Next, let u E H^Q), and define by (19). In addition, let t> (1 ^ be defined 
by 

— / t u ds fc=l»2, ffl) 

hk J (7k 

where <r k is traversed positively and t denotes the unit tangent along a k . Now 
we have: 

Lemma 8.2. The estimate 

|| v (i) _Tu (1) || lv < -^-max(ft,ft , )||u||i i n 


holds. 

Proof. Let cr k denote an edge of r ; - 6 t. The functional 

u) :=(4 1) - t * ' U / J) ) 


is linear, bounded on H 1 (r ; ) by the trace theorem, vanishes on constant u, and 
so has the estimate 

Mu) I < C'( r ;)|u|l,ry • 

By the usual mapping to a reference element [3] it follows that 


C(tj) < 


C 

Q' 


Following the notation of the previous proof, using Ai + A 2 = 1 we have 

(u (1) - Tu (1) ) k < A(M(u)| + M(u)l 

< ^(|u|l.r kl + |u|l,T kj ). 
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and consequently with a modified constant C, 

£ !(,<■> -TVVtf* < (£i)» m «(jr<*>)£ |»|?, TJ 

* i 

which is equivalent to the required result, □ 

Combining lemmas 8.1 and 8.2 now gives: 

Lemma 8.3. Lei it 6 U denote the covolume approximation to u £ H 1 (Q) 
estimated in theorem 6.1 and let Tu denote the approximate tangential compo- 
nents. Then the estimate 

||t> (1 ) — Tu\\w < -1^-max (h, /i')||u|| lin 

holds where i/ 1 ) is defined by (27). 

Proof. It is only necessary to observe that 

l|w (I) - T «\\w < ||v (1) — Tu^Ww + ||T(u — u^)][iy 

and to use lemma 8.2 on the first term and lemma 8.1 and theorem 6.1 on the 
second. □ 

We may now construct the discrete vector field u c „ := («, Tu) G U x U and 
regard it as a covolume approximation to u the solution of the div-curl system 
OH 4 )- The error estimate for ||(u -u^, Tu-vW ) ||^ := ||u- u (1) ||^ + ||Tu- 
v 1 1 vv follows immediately. Defining := (tH), tH) we summarize with 

Theorem 8.1. The covolume approximation H 1 ) satisfies 

ll u cv - u (1) ||w < ^ma x(h, A')ll u l|i.n 

where 0 denotes rnin m |s*n0 m |, 0 m denoting the angles of the triangles in r, A 
denotes the largest absolute interpolation coefficient, and 6 2 denotes the largest 
kite ratio for triangles in the tri angulation. 

Proof. This is simply a question of applying the definition of the norm on the 
left and theorem 6.1 and lemma 8.3. □. 

We can make this bound more explicit for acute angled triangulations. In that 
case, we have A < 1. For any acute angled triangle tj G t, let a[ J) > A ( f > 
A^ > 0 denote the areas of the three triangles formed when the circumcenter 
is joined to the vertices. Thus a\^ f\rj\ are the barycentric coordinates of the 
circumcenter. Obviously, A ^ ■* > |rj|/3. A calculation shows that A^ < 2Ai?\ 
We will also use the area ratio of the triangulation defined by 

| r .| 

r r := max -q (r,, tj adjacent). 
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Note that this ratio is formed for adjacent triangles only. 

Denote the kite areas of tj by I\[^ > > 0 (since h k > 0). Below, 

we will use the fact that for acute triangles we always have < K^\ This 
follows since the maximal property of ensures that it is always at least as 
large as the which are not part of k[*\ But at least one of these two 
must be at least equal to A^\ so the above property does hold. It follows that 


Ki (tj) 


— 

h k tij 2 


< 

A l( T i) + r rl 

r ;l 

< 

2A 2 (rj) + r. 

r 3j4i(r ; ) 

< 

2K 2 (Tj) + r 

, 3 - 2 * 2 ( 73 ) 

< 

8 r T K 2 (Tj). 



This shows that we can take 8 = yJSr T . It would of course be possible to bound 
r r in terms of mesh lengths and 0 but we prefer not to do this. Instead, making 
use of this value for 6 in theorem 8.1 now gives the estimate 

||u«v -u (1) ||fv < ^-^max(/i, ft')||u||i ( n. 

This is valid for acute triangulations. 

Estimating A and 8 is more difficult when there are obtuse triangles. The 
problem is associated with h' k /hk values which approach zero. This will not be 
discussed further here, but section 9 has some related information. 

9 Comments on rectangular meshes. The earlier results have been ob- 
tained under the assumption that all h k > 0. This is equivalent to assuming 
that circumcenters of adjacent triangles are distinct. Now we will briefly summa- 
rize the necessary changes for the contrary case. The previous results continue 
to hold with some slight changes. 

Coincident circumcenters will occur if the triangulation contains a cyclic quadri- 
lateral Q. The distance between the circumcenters of the two triangles making 
Q is h f :== 0. There is no change in the definition of the covolumes, and the 
discrete circulation equations (9) continue to approximate the continuous ones 
(2). Note that the normal component for the diagonal does not appear in the 
circulation equations. The flux equations are also unchanged by the existence 
of Q, but the diagonal component does appear in them. A problem shows up 
in the analysis, since the kite area for the diagonal is zero and consequently 
|| • || w is now merely a seminorm. This inconvenience is easily removed by a 
slight modification of the equations: we simply add the two flux equations asso- 
ciated with Q. The diagonal component does not appear in the resulting sum. 
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In this way, the number of equations is reduced by one and so is the number of 
equations. 

More directly, we can modify the triangulation by deleting the diagonals of all 
cyclic quadrilaterals. Circumcenters are well defined for each of the resulting 
mesh figures. Flux and circulation equations are written in the obvious way. 
Euler’s formula, lemma 2.1 is still true and there will be a single consistency 
condition analogous to (5). The effect of this modification is to remove the dual 
edges which have h f — 0. With these gone, we can prove the key theorem 5.1. 
which requires the existence of ( h f )~ l . Still defining the kite areas by joining 
circumcenters of mesh figures to their vertices, the results in sections 6-7 follow 
essentially verbatim. 

These observations apply particularly to cartesian rectangular meshes. Normal 
to vertical (horizontal) mesh segments we have horizontal (vertical) field compo- 
nents. Flux equations are written for mesh rectangles, and circulation equations 
for the dual mesh rectangles around the mesh points. The two potentials 4> and 
ip are defined respectively at the mesh cell circumcenters and the mesh points 
themselves. The mesh pair is the same as the usual staggered mesh system 
found in fluid flow discretizations. Error estimation of this scheme is given at 
the end of section 6, where second order accuracy is demonstrated. The scheme 
itself is advocated in [2]. 

There are occasions when rectangular and triangular meshes can be usefully 
combined. This is especially true when highly stretched meshes must be used. 
In regions where stretching is necessary, rectangular mesh cells can be used in 
place of triangles. This might be helpful in calculations involving boundary and 
other kinds of layers. 

A situation related to the above occurs when some h' k /hk is small, but positive. 
Although, in a sense, this is inconsequential from the analytical viewpoint, in 
finite precision environments instabilities could result. The three dimensional 
case is worse in this respect, because there is more opportunity for degeneracy. 
In practice, it should be satisfactory to merely set to zero those h! for which 
h' k /hk is below some threshold and proceed as indicated above. This does 
introduce some extra error, but by a suitable choice of the threshold it should 
be within the discretization error. 


10 Mixed boundary conditions. In this section, we will extend the earlier 
results to boundary conditions which are partly normal and partly tangential. 
The technique will be illustrated for the simply connected problem. If required, 
multiply connected cases can be handled by direct extension of the method 
already used in earlier sections. 

We consider in the simply connected polygonal domain Q the system 

div u — p in ft (28) 
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curlu = 

in Q 

(29) 

u • n = / 

on r n 

(30) 

II 

3 

on T t 

(31) 


where T — T n UT t and neither part has zero measure. For illustration, we can 
assume that V is divided into two continuous parts T n and T t each with nonzero 
length. We assume that p E ), w E L 2 (Cl ) and / E i7 1 ^ 2 (T ) and that the 

system has a unique solution uGH 1 ^ ) 

Discretization of (28)-(31) follows section 3, except that (10) is applied only 
to boundary edges on T n . Here and below, it is assumed that nodes of the 
triangulation always separate r n and T t and that, conventionally, the separating 
nodes belong to T t . Boundary edges of r which are in r„ are labelled E\ + 
1, ..., E 2 and boundary edges in T* are labelled + 1, • ••, E. Then in place of 
(10) we have 

u k - 7 - f fds h = Ex + 

hk Ja h 

To complement the boundary conditions, some boundary circulation equations 
are used. The derivation of these equations will be illustrated using Figure 7. 
ABC is part of T t and the boundary covolume r B associated with the node at 



Figure 7: 

B is shown. The orientation conventions are unchanged. Define 

:= mL 9d ‘ 

and 


" BC ' = Wc\I BC 3d ’' 

where AB and BC are the positive directions. The circulation equation for dr f B 
is approximated by 


+ 2 Vab ^ ab + 2 Vbc ^ bc = 



w dx dy 


(32) 
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and a similar equation is used for each boundary covolume associated with a 
node strictly interior to T*. It is always assumed, without explicit mention, 
that h is sufficiently small for to contain at least one strictly interior node. 
This assumption is essential for what follows. If it is not satisfied, then we are 
essentially in the situation of the previous sections as far as the discretization 
is concerned. 

The system of equations which results may be written as 


II 

Ik, 

(33) 

*. 

to* 

1 

•3 

II 

9 

<5 

(34) 


u| r = / 

(35) 

_ 


where g is found from the terms involving the tangential boundary function 
in equations like (32) and C& has no explicit dependence on the lengths of the 
boundary edges. 

The number of nodes interior to T t is E — E^ and the number of edges is therefore 
E — £*2 — 1. In all, the number of equations is 

T 4- Ni + (E ~ E 2 - 1) + (E 2 - Ex) = T + N - 1 

= E 

by lemma 2.1. This time, unlike the earlier situation the number of equations 
and unknowns are equal. 

Now we will prove that the solution of the linear system is unique. For this, 
we will generalize the Helmholtz type theorem 5.3 to cover the new boundary 
conditions. It is necessary first to have results analogous to theorems 4.1 and 
5.1. 

Recalling the definition of the matrix D in section 4, let D* denote the restriction 
of D to {z;; Xi E ^ U FJ, and let B & denote B with columns 1, . . . , E 

deleted. Thus, dim R(Bb ) = E\ + E — E 2 =: E ' . 

Theorem 10.1 Let fi be simply connected and assume that v E R E ' . Then 
there exists <f> E R T such that v = Bb<f> iff D\v — 0. 

Proof. The proof of this follows closely the proof of theorem 4.1 and we shall 
omit most of the details. We have to solve 


Bb<fi = v . 


That 


B\D h $ = 0 

follows from the easily proved relation B ( Drp = 0. Also, rank (D^) = N', where 
N' denotes the number of interior nodes plus nodes on r ( . In addition, 


dimN(Bl) -E' -T. 
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Counting and use of lemma 2.1 shows that N' = E' - T so the result follows. □ 
Now we have that C b u = 0 iff D\H'u = 0. Define S b := (denoting the 

restriction of to W) b y the same s y mbo1 )' 11 follows that CiU = 0 lfF 

u = 5t,<^ for some <p £ R T . 

Next, define 

Uq := {u 6 R E \u\r, = 0} 

W' 0 := {u G f/o;C' 6 u = 0} 

z' := {u G Uq\Fu = 0}. 


Theorem 10.2 Assume that f2 is simply connected. Then Uq has the decompo- 
sition 

Uq = Wq® Z'q 

which is orthogonal relative to [•, •]. 

Proof, tii £ Hj implies [u, iu] = [u.S'i.d 1 ] = [ u . S<j>] = ( Fu , <t > ). Hence, [u, w] = 
0 for all w G W 0 implies that Fu = 0. The rest of the argument is as before. □ 

Existence and uniqueness for the discrete system follows directly, since if u is 
a solution of the homogeneous system we have u G Wq fl Z' Q = 0. Existence 
follows from uniqueness since the coefficient matrix is square. 

It remains to estimate the error of this approximation. The method is similar 
to section 7 with some small modificaions. We define t/ 1 ' just as in (19). Cor- 
responding to u we define a very slightly different function, still denoted by 
as follows: 

u (2) := _L [ u .nds *=1,2,. ...Ei k = E 2 + l,...,E 

K J*'„ 

4 2) : = k=E^ + l,...,E 2 . 

The edges <r* k = E 2 + 1, • ••, E are in T t ; each of the coedges a' k extends from 
a k ’s midpoint to the circumcenter of the triangle containing a k , exemplified by 
PQ and TU in Figure 7. As before, define e (<) -u- u (,) * = 1,2. Then 

Fe w = 0 e Z' Q 
C b eW = 0 e< 2) G W' 0 

and [cO), €«] = 0 by theorem 10.2. Following the method of section 6, we 
obtain, corresponding to (20) 

||u-u (1) ||vv < ||u (2) - u^Ww- 
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The right side can be estimated as before, except that now there will be a 
contribution to the norm from the boundary segments along IV The earlier 
estimation technique works here too, and is especially simple if each coedge 
associated with lies wholly in This will certainly be the case, for example, if 
the boundary triangles are acute angled, since then the coedges will lie inside the 
triangles themselves. That is the simplest case. We summarize this discussion 
with: 

Theorem 10.3 For the equations ( 28)- (31) and their discretization ( 33 )-(S5 ), 
we have the estimate 


||« - < Cmax(/»,/i')|«| lin . 


So far, we have considered the purely normal boundary condition and the mixed 
boundary condition. For the purely tangential case the simplest approach seems 
to be to use tangential field components in place of the normal ones used in 
sections 2-6. The circulation equations are formed using the triangle boundaries 
and the flux equations are formed using the covolume boundaries. It is apparent 
that this method gives a natural dual to the earlier one. We will not go into the 
details here, but an almost verbatim development of the theory can be carried 
out. This includes most of the results of sections 2-6. The main difference is to 
exchange the roles of <p and ip. This duality is strongly reminiscent of complex 
variable theory. In fact, a similar connection was developed in [5]. The results 
given above provide useful tools for the analysis of that scheme. 
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